function [err, EPSCresult]=EPSC_MeanVar()
% Analyze the Mean and Variance of a set of EPSPs
% Derived from EPSP_COV
% Expects data to be a series of EPSCs in a train, 10 msec apart.
% Oleskovich et al., 1999
% Paul B. Manis, Ph.D.
% pmanis@med.unc.edu
% Initial versions: 4/2002, 5/2002.
% 6/2/2002.
% 12/1/02
%

try
   sf = getmainselection;
   if(sf > 0) 
      pflag = getplotflag;
      QueMessage('EPSC_MeanVar analysis', 1); % clear the que
      for i = 1:length(sf)
         EPSC_MeanVar2(sf(i), pflag);
      end;
   end;
catch
   watchoff;
   QueMessage('Error in EPSC_MeanVar Analysis routine', 1);
end;

%-----------------------------------------------------------

function [err, EPSCresult]=EPSC_MeanVar2(sf, plot_flag, pflag)

global VOLTAGE CURRENT DFILE
global CONTROL
% first generate empty output arrays

err = 0;
EPSCresult=[];

[DFILE, err] = analysis_setup(DFILE, sf); % do default setup.
if(err ~= 0)
   return;
end;

%-----------------------------Begin analysis--------------------------------
QueMessage('EPSC_MeanVar - Starting');
dat = [];
time = [];

% abstract general information
protocol=deblank(lower(CONTROL(sf).protocol));
rate = DFILE.rate.*DFILE.nr_channel/1000;
[records,pts]=size(CURRENT);

[RL, err] = record_parse(CONTROL(sf).reclist);
frec=min(RL);
lrec=max(RL);

QueMessage('EPSC_MeanVar - time base');
% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
k = find(DFILE.rate < 1);
DFILE.rate(k) = 50;
RATES = (DFILE.rate .* DFILE.nr_channel)/ 1000; % array of sampling rates, convert to msec

% always compute the time of each stimulus, in seconds from the start of the data
QueMessage('EPSC_MeanVar - time base 2');
if(DFILE.mode >= 5)
   wz=DFILE.ztime;
   w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
   if(~isempty(w))
      wz(w+1:end)=wz(w+1:end)+24*60*60;
   end;
   zt = (wz-wz(1))/(60);
   cond_baseline = 5; % 5 min baseline
else
   zt = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
   cond_baseline = 3;
end;
TM=zt;

QueMessage('EPSC_MeanVar - time base 3');
if(DFILE.mode >= 5)
   wz=DFILE.ztime;
   w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
   if(~isempty(w))
      wz(w+1:end)=wz(w+1:end)+24*60*60;
   end;
   ZT = (wz-wz(1))/(60);
   cond_baseline = 5; % 5 min baseline
elseif length(DFILE.ztime) == length(RL)
   ZT = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
   cond_baseline = 3;
else
   ZT = ones(length(RL), 1)*DFILE.cycle;
   cond_baseline = 5;
end;

% access ZT with the record number (1..n) to get the corresponding ztime
% have to get the records first

% Now, get the times when the valves switched (if any...)
% and generate periods with valve 1, 2, 3 or 4.
QueMessage('EPSC_MeanVar - Valves...');
p=datac('getnote'); % read the current notefile information.
t_sw_valve=[];
n_valve=[];
TL=[];
VL=[];
if(~isempty(p) & length([p.proto]) > 0) % there should be some, but if not, don;'t do much
   % first set of arrays are immediate representations.
   sw_valve=[1 [p(find(diff([p.valve])~=0)+1).frec]];	% valve switch list (records)
   if(length(sw_valve) > 1)
      n_valve=[1 p(find(diff([p.valve])~=0)+1).valve]; % which valve...
      t_sw_valve=ZT(sw_valve); % don't forget offset from start of data..
      % now make long time arrays to match the other arrays.
      for i = frec:lrec
         TL(i)=ZT(i);
         for j=1:length(sw_valve)
            if(i >= sw_valve(j))
               VL(i)=n_valve(j);
            end;
         end;
      end;
   end;
end;

% Get analysis windows 
QueMessage('EPSC_MeanVar - analysis windows');
stim_list=CONTROL(sf).stim_time; % get the array
psp_time=CONTROL(sf).psp_time; % get psp definition array
if(ischar(psp_time))
   psp_time = str2num(psp_time);
end;
if(ischar(stim_list))
   stim_list= str2num(stim_list);
end;
%l = length(psp_time);
%if(l > 4)
%	psp_time = reshape(psp_time, l/2, 2)';
%end;
p2=size(psp_time);
npsc = length(stim_list); % get number of pscs to examine

trmp = floor(number_arg(stim_list(1))./RATES - 0.5); % for RMP/ihold determination
for i = 1:npsc % compute time windows for each psc in the train
   t0(i,:) = floor(number_arg(psp_time((i-1)*2+1))./RATES);
   t1(i,:) = floor(number_arg(psp_time((i-1)*2+2))./RATES);
end;

%-----------------------------Begin analysis in earnest --------------------------------
QueMessage('EPSC_MeanVar - Measuring Ih and rmp');

% measure holding current and "rmp"
% also measure Rin with s1 pulse after ho (ts2)
% and tau by fitting to s1 pulse.
for i = 1:records
   hold_cur(i) = mean(CURRENT(i,1:trmp(i))); % save in array for later usage
   RMP(i) = mean(VOLTAGE(i,1:trmp(i)));
end;
CONTROL(sf).iHold = mean(hold_cur);
CONTROL(sf).Rmp = mean(RMP);
[CURRENT] = FP_artsupp(CURRENT, DFILE, sf);
% Now smooth the Current out a bit
QueMessage('EPSC_MeanVar - Smoothing');
for i = 1:records
   fsamp = 1000/RATES(i); % get sampling frequency
   fco = 10000;		% cutoff frequency in Hz
   wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
   if(wco < 1) % if wco is > 1 then this is not a filter!
      [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
      ismo(i,:) = DigitalFilt(b, a, CURRENT(i,:)')'; % filter all the traces...
   else
      ismo(i,:) = CURRENT(i,:);
   end
end

QueMessage('EPSC_MeanVar - Finding EPSC peak currents');
for r = 1:records
   ibase(r)=mean(CURRENT(r,1:trmp(r)));
end;
% find EPSC amplitudes for each entry in the train

% get psc amplitudes for every pulse and record.
for r = 1:records
   ismo(r,:) = ismo(r,:) - ibase(r);
   for i = 1:npsc
      I_psc(i, r) = min(ismo(r,t0(i,r):t1(i,r)))';
   end;
   zvar(r) = ismo(r, 5); % get the variance from a single point measurement before the stimulus
end;

% reduce to a list without failures - must be true for all pulses in consequective records.
r_ok = [];
failthresh = 200; % in pA - this is biggern than minis in normal conditions.
for r = 1:records-1
   if(all(abs(I_psc(:, r)) > failthresh) & all(abs(I_psc(:, r+1)) > failthresh))
      r_ok = [r_ok r];
   end;
end;

mzvar = var(zvar(r_ok)); % get the baseline noise variance - use to subtract.

% compute mean current for each stimulus pulse
for i = 1:npsc
   I_m(i) = mean(I_psc(i,r_ok));
   I_var(i) = var(I_psc(i,r_ok))-mzvar;
end
ipscmin = min(min(I_psc));

Cvq = 0.32; % quantal variability; is CV of mEPSCs - use to correct per Scheuss, eq. 7.

% fit the parabola to the 
EPSC_MeanVar.mzvar = mzvar;
EPSC_MeanVar.Im = I_m;
EPSC_MeanVar.Ivar = I_var;
EPSC_MeanVar.npsc = npsc;

CONTROL(sf).EPSC_MeanVar = EPSC_MeanVar; % copy the structure over

QueMessage('EPSC_MeanVar - analysis complete');


%-----------------------------Prepare for plotting--------------------------------
% for plotting, do baselines/stddev.
%----------------------------- plot if figure is set--------------------------------
if(plot_flag >= 0)
   h = findobj('Tag', 'EPSC_MeanVar'); % check for pre-existing window
   if(isempty(h)) % if none, make one
      h = figure('Tag', 'EPSC_MeanVar', 'Name', 'EPSP Time Course Analysis', 'NumberTitle', 'off');
   end
   figure(h); % otherwise, select it
   clf; % always clear the window...
   fsize = 7;
   msize = 3;
   tmax = max(time);    
   
   % plot Max voltage for EPSP
   subplot('Position',[0.53,0.75,0.20,0.20]);
   sym = {'ks', 'ko', 'k^', 'k+', 'ks'};
   for i = 1:npsc
      plot(r_ok, -I_psc(i,r_ok)/1000, sym{i}, 'Markersize', 2, 'Markerfacecolor', 'k') % data in black
   end;
   
   % plot the histogram of amplitudes (are there multiple clear peaks?)
   NBINS = floor(length(r_ok/10));
   [nh, xh] = hist(-I_psc(i,r_ok), NBINS); % 64 bins...
   subplot('Position', [0.78 0.75 0.20 0.20]);
   cla;
   bar(xh, nh);
   set(gca, 'fontsize', fsize);
   xlabel('Amplitude (pA)');
   ylabel('# events');
   % datac_setcrosshair(gca, 'Tr2_hist', 'pA', 'N', [0.055 0.26 0.1 0.07]);
   set(gca, 'FontSize', fsize);
   ylabel('EPSP Amplitude (nA)', 'FontSize', fsize);
   set(gca, 'XTickLabelMode', 'Manual');
   %grid;
   u=get(gca, 'YLim');
   
   % plot the EPSPs...
   subplot('Position', [0.07, 0.55, 0.45, 0.4]);
   plot(time(r_ok,(t0(1,1):t1(1:1)))', ismo(r_ok,(t0(1,1):t1(1:1)))'/1000);
   set(gca, 'XLim', [6 8]);
   set(gca, 'YLim', [ipscmin*1.1/1000 0.5]);
   xlabel('ms', 'FontSize', fsize);
   ylabel('nA', 'FontSize', fsize);
   set(gca, 'FontSize', fsize);
   
   set(gca, 'FontSize', fsize);
   ylabel('EPSP Amplitude (nA)', 'FontSize', fsize);
   set(gca, 'XTickLabelMode', 'Manual');
   %grid;
   u=get(gca, 'YLim');
   
   % plot the difference between the EPSPs and the mean EPSP...
   ism = mean(ismo);
   j=1;
   for i = r_ok
      ivm(j,:) = ismo(i,(t0(1,1):t1(1:1))); % 
      ivx(j,:) = ism(1,(t0(1,1):t1(1:1)))-ivm(j,:);
      tvx(j,:) = time(i, (t0(1,1):t1(1:1)));
      j = j + 1;
   end;
   subplot('Position', [0.07, 0.1, 0.45, 0.4]);
   plot(tvx', ivx');
   set(gca, 'XLim', [6 8]);
   % set(gca, 'YLim', [ipscmin*1.1/1000 0.5]);
   xlabel('ms', 'FontSize', fsize);
   ylabel('nA', 'FontSize', fsize);
   set(gca, 'FontSize', fsize);
   
   % plot textual information abstracted from analysis...
   %subplot('Position',[0.55,0.1,0.4,0.45])
   	axis([0,1,0,1])
   axis('off')
   text(0,0.9,sprintf('%-12s  R[%d:%d]  Protocol:  %-15s',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol), 'Fontsize', 8, 'Interpreter', 'none');
   text(0, 0.8, sprintf('Comment:  %s', DFILE.comment), 'Fontsize', 8, 'Interpreter', 'none');
   text(0,0.7,sprintf('Sol:%-12s   Gain:%4.1f  LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 8);
   text(0, 0.6, sprintf('mean:  %7.3f   variance:  %8.3f  pA  noise variance: %8.5f pA',...
      I_m(1)/1000, I_var(1)/1000^2, mzvar(1)/1000^2), 'Fontsize', 8);
   
   orient landscape
   drawnow
   % control printing and closing of window for automatic runs
   % f = 1 creates plot and leaves it up
   % f = 2 creates plot but closes it when done
   % f = 3 creates plot and prints it and then closes it
   if (plot_flag > 0)
      print -dljet3;
   end
   if plot_flag == 2  
      close
   end
end;

